Shifts in bird and butterfly communities of Tamhini Ghat

Introduction

This file was last update on 2021-06-30

#Required libraries

library(tidyverse)
library(broom)
library(lubridate)
library(viridis)
library(RColorBrewer)

#Data import

abn <- read.csv("./Data/Abundance_280820.csv")
transects <- read.csv("./Data/Transects.csv")

# adding relevant variables

abn <- abn%>%
  mutate(Year = year(dmy(Date)))%>% #adding year of sampling
  mutate(Study = case_when(Year %in% c(1998, 1999, 2000, 2001) ~ "1998 - 2000", #previous study
                           Year %in% c(2016, 2017) ~ "2016 - 2017"))%>% #current study
  mutate(Transect = as.character(Transect),
         Transect = ifelse(Transect == "T6A" | Transect == "T6B", "T6", Transect),
         Transect = as.factor(Transect))

transects <- transects%>%
  mutate(Transect = as.factor(Transect))

# cleaning abundance
abn <- abn%>%
  filter(!is.na(Abundance) &
         Specific.Epithet != "" & #removing samples with no species
         Specific.Epithet != " " &
          Specific.Epithet != "Unidentified" & #removing unidentified species
         substr(Specific.Epithet, 1, 1) != "?")%>% #83 obs removed
  droplevels()

theme_set(theme_bw()+
            theme(legend.position = "top",
                  text = element_text(size = 20)
                  )
          )#setting ggplot theme

# function for t-test reporting

tstats <- function(x){
  x%>%
    rename('1998-2001' = estimate1,
         '2016-2017' = estimate2)%>%
    dplyr::select(`1998-2001`, `2016-2017`, parameter, statistic, p.value)
}

Sampling effort

No. of visits to each transect/site.

abn%>%
  group_by(Study, Taxa, Transect)%>%
  summarise(n = length(unique(Date)))%>%
  spread(Transect, n)
Study Taxa T1A T1B T2 T3 T4 T5 T6
1998 - 2000 Birds 19 18 7 10 4 1 1
1998 - 2000 Butterflies 14 7 4 9 4 1 1
2016 - 2017 Birds 23 23 22 19 16 22 17
2016 - 2017 Butterflies 23 20 20 22 18 21 15

No. of individuals sampled

abn%>%
  group_by(Study, Taxa, Transect)%>%
  summarise(n = sum(Abundance))%>%
  spread(Transect, n)
Study Taxa T1A T1B T2 T3 T4 T5 T6
1998 - 2000 Birds 255 165 103 213 236 7 28
1998 - 2000 Butterflies 93 59 12 180 157 10 4
2016 - 2017 Birds 480 393 358 289 165 203 133
2016 - 2017 Butterflies 432 188 207 652 187 224 124

Uneveness in sampling effort across sites in studies may bias estimates of change.

Rarefying across sites

library(vegan)

# input for vegan

veg_birds <- abn%>%
  filter(Taxa == "Birds")%>%
  group_by(Study, Transect, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = paste(Study, "_", Transect))%>%
  group_by(sample)%>%
  #mutate(N = sum(n, na.rm = T))%>%
  #filter(N > 0)%>%#removing samples where no species were detected
  dplyr::select(sample, Specific.Epithet, n)%>%
  ungroup()%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

veg_butterfly <- abn%>%
  filter(Taxa == "Butterflies")%>%
  group_by(Study, Transect, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = paste(Study, "_", Transect))%>%
  group_by(sample)%>%
  #mutate(N = sum(n, na.rm = T))%>%
  #filter(N > 0)%>%#removing samples where no species were detected
  dplyr::select(sample, Specific.Epithet, n)%>%
  ungroup()%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

# sample info

samp <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  dplyr::select(Study, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))

## Species accumulation curves

rarebutts <- rarecurve(veg_butterfly) ## butterflies

# Coerce data into "long" form.
## taken from https://stat.ethz.ch/pipermail/r-sig-ecology/2018-December/005867.html

names(rarebutts) <- samp$sample

protox <- mapply(
  FUN = function(x, y) {
    mydf <- as.data.frame(x)
    colnames(mydf) <- "value"
    mydf$sites <- y
    mydf$subsample <- attr(x, "Subsample")
    mydf}, 
  x = rarebutts, y = as.list(names(rarebutts)), SIMPLIFY = FALSE)

rarebutts.df <- do.call(rbind, protox)

rownames(rarebutts.df) <- NULL  # pretty

rarebutts.df$Taxa <- "Butterflies"

## repeating for birds

rarebirds <- rarecurve(veg_birds)

# Coerce data into "long" form.
## taken from https://stat.ethz.ch/pipermail/r-sig-ecology/2018-December/005867.html

names(rarebirds) <- samp$sample

protox <- mapply(
  FUN = function(x, y) {
    mydf <- as.data.frame(x)
    colnames(mydf) <- "value"
    mydf$sites <- y
    mydf$subsample <- attr(x, "Subsample")
    mydf}, 
  x = rarebirds, y = as.list(names(rarebirds)), SIMPLIFY = FALSE)

rarebirds.df <- do.call(rbind, protox)

rownames(rarebirds.df) <- NULL  # pretty

rarebirds.df$Taxa <- "Birds"

## joining rarefaction values

rarefied <- full_join(rarebirds.df, rarebutts.df)

## adding sample info

rarefied <- left_join(rarefied, samp, by = c("sites" = "sample"))
# Plot.
ggplot(rarefied, aes(x = subsample, y = value, col = Study))+
  geom_line(size = 1)+
  labs(x = "Sample Size", y = "Species")+
  facet_grid(Taxa ~ Transect)

ggsave(last_plot(), filename = "./Figures and Tables/figureA.png", height = 8, width = 16, dpi = 300)  

Change in diversity across years

Species richness

Across the two studies

abn%>%
  group_by(Taxa, Study)%>%
  summarise(richness = length(unique(Specific.Epithet)),
            N = sum(Abundance, na.rm = T))
Taxa Study richness N
Birds 1998 - 2000 70 1007
Birds 2016 - 2017 105 2021
Butterflies 1998 - 2000 45 515
Butterflies 2016 - 2017 66 2014

More species of both taxa were encountered in current study compared to the previous study. Difference in sampling effort is apparent accross two studies.

Across sites in the two studies

No. of species encountered across sites sampled in the two studies.

abn%>%
  left_join(transects, by = c("Transect"))%>%
  group_by(Taxa, Study, Transect)%>%
  summarise(richness = length(unique(Specific.Epithet)))%>%
  spread(Transect, richness)
Taxa Study T1A T1B T2 T3 T4 T5 T6
Birds 1998 - 2000 47 37 28 36 29 4 11
Birds 2016 - 2017 66 59 53 47 31 38 28
Butterflies 1998 - 2000 27 22 8 37 24 5 3
Butterflies 2016 - 2017 41 32 36 44 31 39 27

Fewer species were encountered in sites with lower effort in the previous study, i.e., T5 and T6.

Shannon diversity

# Function to caluculate Shannon diversity

shanon <- function(c){
  ## c is a vector of species abundances in a site 
  ## Spade cannont compute n < 20
  
  if(sum(c) < 20){
    
    H <- "Sample size too small"
    
  }else{
    
    require(SpadeR)
    require(tidyr)
    
    div <- Diversity(c, datatype = "abundance")
    
    H <- data.frame(div$Shannon_diversity)%>%
      rownames_to_column("Method")
    
  }
  
  return(H)
  
}

# Calculating diversity in each site over years

div_year <- abn%>%
  group_by(Study, Taxa, Transect,  Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  spread(Specific.Epithet, n, fill = c(0))%>%
  group_by(Study, Taxa, Transect)%>%
  nest()%>%
  mutate(H = map(data, ~shanon(.)))%>%
  dplyr::select(-data)%>%
  unnest(H)%>%
  filter(grepl("MLE", Method) | is.na(Method))

# Summary

div_year%>%
  group_by(Taxa, Study)%>%
  skimr::skim(Estimate)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

Variable type: numeric

Taxa Study mean sd
Birds 1998 - 2000 17.31 6.47
Birds 2016 - 2017 24.88 5.74
Butterflies 1998 - 2000 15.87 3.42
Butterflies 2016 - 2017 20.40 3.24

Linear model testing change in diversity across years:

div_year%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~lm(Estimate ~ Study, data = .)),
         sumry = map(test, broom::tidy),
         stat = map(test, broom::glance))%>%
  dplyr::select(sumry, stat)%>%
  unnest()%>%
  select(Taxa:r.squared)
## Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'select' for signature '"grouped_df"'

Plotting diversity across studies:

div_year%>%
  ggplot(aes(Taxa, Estimate, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Effective number of species")+
  scale_fill_brewer(palette = "Greys")

ggsave(last_plot(), filename = "./Figures and Tables/figure2.png", height = 6, width = 8, dpi = 300)  

The first order diversity of both taxa has increased across the two surveys. SIgnificantly for birds but not for butterflies. Likely due to low sampling effort.

Change in composition of species accross years

Turn over of each taxa across studies:

# sample info

samp_birds <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  filter(Taxa == "Birds")%>%
  dplyr::select(Study, Taxa, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))

samp_butterfly <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  filter(Taxa == "Butterflies")%>%
  dplyr::select(Study, Taxa, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))

# Similarity matrix by studies

sim_birds <- abn%>%
  filter(Taxa == "Birds")%>%
  group_by(Study, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = Study)%>%
  ungroup()%>%
  dplyr::select(sample, Specific.Epithet, n)%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

sim_butterfly <- abn%>%
  filter(Taxa == "Butterflies")%>%
  group_by(Study, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = Study)%>%
  ungroup()%>%
  dplyr::select(sample, Specific.Epithet, n)%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

# bray - curtis

data_frame(
  Taxa = c("Birds", "Butterflies"),
  Dissimilarity = c(as.numeric(vegdist(sim_birds, method = "bray")),
                    as.numeric(vegdist(sim_butterfly, method = "bray")))
)
# ordination

NMDS_birds <- metaMDS(veg_birds, distance = "bray", k = 2, model = "local")

NMDS_butterfly <- metaMDS(veg_butterfly, distance = "bray", trymax = 1000, k = 2, model = "local")


# Extracting scores and adding habitat vars

ord_birds <- data.frame(scores(NMDS_birds))%>%
  rownames_to_column("sample")%>%
  left_join(samp_birds, "sample")

ord_butterfly <- data.frame(scores(NMDS_butterfly))%>%
  rownames_to_column("sample")%>%
  left_join(samp_butterfly, "sample")

ord <- bind_rows(ord_birds, ord_butterfly)

Testing change in composition of birds :

adonis(veg_birds ~ Study, data = samp_birds)
## 
## Call:
## adonis(formula = veg_birds ~ Study, data = samp_birds) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## Study      1   0.77547 0.77547  3.9632 0.24827  0.002 **
## Residuals 12   2.34798 0.19567         0.75173          
## Total     13   3.12345                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Testing change in composition of butterflies:

adonis(veg_butterfly ~ Study, data = samp_butterfly)
## 
## Call:
## adonis(formula = veg_butterfly ~ Study, data = samp_butterfly) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Study      1    0.9191 0.91908   4.054 0.25253  0.001 ***
## Residuals 12    2.7205 0.22671         0.74747           
## Total     13    3.6396                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

NMDS ordination plot:

# plotting

ord%>%
  ggplot(aes(NMDS1, NMDS2, col = Study, fill = Study))+
  stat_ellipse(level = 0.95, geom = "polygon", linetype = "dashed", size = 1, alpha = 0.3)+
  geom_point(aes(shape = Habitat), size = 3)+
  guides(shape = guide_legend(nrow = 2, byrow = T))+
  scale_color_brewer(palette = "Set1")+
  scale_fill_brewer(palette = "Set1")+
  facet_wrap(~Taxa)+
  theme(legend.box = "vertical")

ggsave(last_plot(), filename = "./Figures and Tables/figure3.png", height = 6, width = 8, dpi = 300)

Both bird and butterfly communities compositions are significantly different compared to the previous study.

Change in niche width of community

Trophic Guilds

Diversity across trophic guilds:

# compiling trophic guild data

host <- read.csv("./Data/Butterflies_host plant_280820.csv") #butterfly host plant data

diet <- read.csv("./Data/Birds_diet.csv") #bird diet data

trophic <- host%>%
  full_join(diet)%>%
  dplyr::select(Specific.Epithet, Guild)%>%
  distinct()

# Calculating diversity of guilds

guilds <- abn%>%
  left_join(trophic, by = c("Specific.Epithet"))%>%
  filter(!Guild == "")%>%
  group_by(Study, Taxa, Transect, Guild,  Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  spread(Specific.Epithet, n, fill = c(0))%>%
  group_by(Study, Taxa, Transect, Guild)%>%
  nest()%>%
  mutate(H = map(data, ~shanon(.)))%>%
  dplyr::select(-data)%>%
  unnest(H)%>%
  filter(grepl("MLE", Method) | is.na(Method))

# summary

guilds%>%
  group_by(Taxa, Guild, Study)%>%
  skimr::skim(Estimate)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Guild, Study, mean, sd)

Variable type: numeric

Taxa Guild Study mean sd
Birds Carnivore 1998 - 2000 3.43 0.31
Birds Carnivore 2016 - 2017 6.26 0.15
Birds Frugivore 1998 - 2000 4.11 0.92
Birds Frugivore 2016 - 2017 4.21 0.60
Birds Grainivore 1998 - 2000 3.86 0.75
Birds Grainivore 2016 - 2017 3.84 0.49
Birds Insectivore 1998 - 2000 8.77 1.95
Birds Insectivore 2016 - 2017 13.58 4.82
Birds Omnivore 1998 - 2000 4.82 0.92
Birds Omnivore 2016 - 2017 5.88 1.15
Butterflies Generalist 1998 - 2000 3.16 1.48
Butterflies Generalist 2016 - 2017 5.95 1.87
Butterflies Grass Specialist 1998 - 2000 1.80 NA
Butterflies Grass Specialist 2016 - 2017 2.71 0.71
Butterflies Herb Specialist 1998 - 2000 3.58 NA
Butterflies Herb Specialist 2016 - 2017 4.09 0.86
Butterflies Liana Specialist 1998 - 2000 NaN NA
Butterflies Shrub Specialist 1998 - 2000 4.36 0.77
Butterflies Shrub Specialist 2016 - 2017 4.18 2.10
Butterflies Tree Specialist 1998 - 2000 5.92 3.07
Butterflies Tree Specialist 2016 - 2017 6.88 1.63

Linear model testing change in guild diversity:

guilds%>%
  filter(Guild != "Liana Specialist")%>% #removing due to lack of obs
  complete(nesting(Study, Taxa), Transect, Guild, fill = list(H = 0))%>%
  group_by(Taxa, Guild)%>%
  nest()%>%
  mutate(test = map(data, ~lm(Estimate ~ Study, data = .)),
         sumry = map(test, broom::tidy),
         stat = map(test, broom::glance))%>%
  dplyr::select(sumry, stat)%>%
  unnest()%>%
  dplyr::select(Taxa:r.squared)
Taxa Guild term estimate std.error statistic p.value r.squared
Birds Carnivore (Intercept) 3.4273333 0.0482458 71.0390510 0.0000000 0.9745403
Birds Carnivore Study2016 - 2017 2.8301667 0.0849454 33.3174666 0.0000000 0.9745403
Birds Frugivore (Intercept) 4.1146000 0.1159192 35.4954216 0.0000000 0.0049246
Birds Frugivore Study2016 - 2017 0.0966857 0.1517738 0.6370381 0.5258746 0.0049246
Birds Grainivore (Intercept) 3.8566667 0.1170495 32.9490144 0.0000000 0.0002457
Birds Grainivore Study2016 - 2017 -0.0166667 0.1850716 -0.0900553 0.9287877 0.0002457
Birds Insectivore (Intercept) 8.7660000 0.6795468 12.8997745 0.0000000 0.2877674
Birds Insectivore Study2016 - 2017 4.8127143 0.8628491 5.5777010 0.0000003 0.2877674
Birds Omnivore (Intercept) 4.8210000 0.1665640 28.9438346 0.0000000 0.2223562
Birds Omnivore Study2016 - 2017 1.0560000 0.2180834 4.8421836 0.0000060 0.2223562
Butterflies Generalist (Intercept) 3.1625000 0.5157184 6.1322230 0.0000001 0.3200482
Butterflies Generalist Study2016 - 2017 2.7838333 0.5738384 4.8512495 0.0000124 0.3200482
Butterflies Grass Specialist (Intercept) 1.7980000 0.2732411 6.5802684 0.0000001 0.2042940
Butterflies Grass Specialist Study2016 - 2017 0.9124000 0.2921071 3.1235123 0.0034125 0.2042940
Butterflies Herb Specialist (Intercept) 3.5810000 0.3907790 9.1637464 0.0000000 0.0300820
Butterflies Herb Specialist Study2016 - 2017 0.5111429 0.4064163 1.2576829 0.2142349 0.0300820
Butterflies Shrub Specialist (Intercept) 4.3575000 0.4412280 9.8758466 0.0000000 0.0024700
Butterflies Shrub Specialist Study2016 - 2017 -0.1781000 0.5220680 -0.3411433 0.7345161 0.0024700
Butterflies Tree Specialist (Intercept) 5.9220000 0.4125034 14.3562448 0.0000000 0.0525717
Butterflies Tree Specialist Study2016 - 2017 0.9577143 0.4930359 1.9424840 0.0562211 0.0525717

Plotting guild diversity across years:

guilds%>%
  ggplot(aes(reorder(Guild, Estimate), Estimate, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5, position = position_dodge(preserve = "single"))+
  labs(y = "Effective number of species", x = "Trophic Guild")+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  scale_fill_brewer(palette = "Greys")+
  facet_wrap(~Taxa, scales = "free_x", ncol = 1)

ggsave(last_plot(), filename = "./Figures and Tables/figure4.png", height = 9, width = 8, dpi = 300)  

Carnivorous, Grainivorous birds have reduced in the comunity, compared to insectivorous birds which have increased. Similarly, Shrub specialist and tree specialist butterflies have reduced where grass specialists and herb specialist have increased. Suprsingly, generalist species of both taxa appear to be unchanged.

Change in niche specialisation

-Can we incorporate abunndance into community specialisation index?

Trophic niche

Trophic specialisation based on Juliard et al. 2006

# 

butterflies_TSI <- host%>%
  mutate(H = length(unique(Hostplant.Family)))%>%
  group_by(Specific.Epithet)%>%
  mutate(h = length(unique(Hostplant.Family)))%>%
  ungroup()%>%
  mutate(TSI = sqrt((H/h)-1))

birds_TSI <- diet%>%
  mutate(H = length(unique(Prey)))%>%
  group_by(Specific.Epithet)%>%
  mutate(h = length(unique(Prey)))%>%
  ungroup()%>%
  mutate(TSI = sqrt((H/h)-1))

CSI.T <- abn%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
  group_by(Study, Taxa, Transect)%>%
  distinct()%>%
  left_join(butterflies_TSI, by = c("Specific.Epithet"))%>%
  left_join(birds_TSI, by = c("Specific.Epithet"))%>%
  mutate(TSI = ifelse(is.na(TSI.x), TSI.y, TSI.x))%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet, TSI)%>%
  summarise(CSI.T = mean(TSI, na.rm = T))

# summary

CSI.T%>%
  group_by(Taxa, Study)%>%
  skimr::skim(CSI.T)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

Variable type: numeric

Taxa Study mean sd
Birds 1998 - 2000 2.59 0.16
Birds 2016 - 2017 2.53 0.09
Butterflies 1998 - 2000 5.34 0.39
Butterflies 2016 - 2017 5.07 0.24

Linear model testing difference is trophic specialization:

CSI.T%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~lm(CSI.T ~ Study, data = .)),
         sumry = map(test, broom::tidy),
         stat = map(test, broom::glance))%>%
  dplyr::select(sumry, stat)%>%
  unnest()%>%
  dplyr::select(Taxa:r.squared)
Taxa term estimate std.error statistic p.value r.squared
Birds (Intercept) 2.5880088 0.0493181 52.4758739 0.0000000 0.051883
Birds Study2016 - 2017 -0.0565188 0.0697463 -0.8103492 0.4335164 0.051883
Butterflies (Intercept) 5.3385626 0.1229059 43.4361931 0.0000000 0.170905
Butterflies Study2016 - 2017 -0.2733716 0.1738151 -1.5727722 0.1417534 0.170905
CSI.T%>%
  ggplot(aes(Taxa, CSI.T, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Community specialisation index \n (Trophic)")+
  scale_fill_brewer(palette = "Greys")+
  scale_y_sqrt()

Niche width of the community did not change over the two studies.

Habitat Niche

HSI <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  group_by(Study, Taxa, Habitat, Specific.Epithet)%>%
  # calculating abundance of each species in each habitat type
  summarise(n = sum(Abundance, na.rm = T))%>%
  group_by(Study, Taxa)%>%
  mutate(N = sum(n), # total number of indviduals encountered 
         rel.prop = n/N)%>% # relative proportion of species
  group_by(Study, Specific.Epithet)%>%
  summarise(HSI = sd(rel.prop)/mean(rel.prop))%>%
  replace_na(replace = list(HSI = 0))

CSI.H <- abn%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
  group_by(Study, Taxa, Transect)%>%
  distinct()%>%
  left_join(HSI, by = c("Study", "Specific.Epithet"))%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet, HSI)%>%
  summarise(CSI.H = mean(HSI, na.rm = T))

# summary

CSI.H%>%
  group_by(Taxa, Study)%>%
  skimr::skim(CSI.H)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

Variable type: numeric

Taxa Study mean sd
Birds 1998 - 2000 0.54 0.10
Birds 2016 - 2017 0.49 0.02
Butterflies 1998 - 2000 0.61 0.09
Butterflies 2016 - 2017 0.57 0.03

Linear model testing for habitat specialisation:

CSI.H%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~lm(CSI.H ~ Study, data = .)),
         sumry = map(test, broom::tidy),
         stat = map(test, broom::glance))%>%
  dplyr::select(sumry, stat)%>%
  unnest()%>%
  dplyr::select(Taxa:r.squared)
Taxa term estimate std.error statistic p.value r.squared
Birds (Intercept) 0.5363977 0.0278655 19.249520 0.0000000 0.0869142
Birds Study2016 - 2017 -0.0421175 0.0394078 -1.068760 0.3062059 0.0869142
Butterflies (Intercept) 0.6110527 0.0244089 25.033991 0.0000000 0.1087676
Butterflies Study2016 - 2017 -0.0417742 0.0345194 -1.210166 0.2495141 0.1087676
CSI.H%>%
  ggplot(aes(Taxa, CSI.H, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Community specialisation index \n (Habitat)")+
  scale_fill_brewer(palette = "Greys")

Habitat niche of the community did not change over the study periods.

Seasonal changes in relative abundance

season <- abn%>%
  mutate(Month = month(dmy(Date), label = T))%>%
  group_by(Study, Month, Taxa)%>%
  summarise(n = sum(Abundance, na.rm=T))%>%
  group_by(Study, Taxa)%>%
  mutate(N = sum(n),
         rel.prop = n/N)

season%>%
  ggplot(aes(Month, rel.prop, col = Taxa, shape = Taxa))+
  geom_point(size = 3)+
  geom_path(aes(group = Taxa), size = 1, linetype = "dashed")+
  scale_color_brewer(palette = "Set1")+
  facet_wrap(~Study, ncol = 1)

Site Map

library(raster)
library(sf)

# Raw point data

points <- read.csv("./Data/sites.csv")%>%
  drop_na()%>%
  dplyr::select(-Location)

# Converting to SF points

points_sf <- st_as_sf(points, coords = c('Lon', 'Lat'))

st_crs(points_sf) <- 4326

# Making transect lines

lines <- points_sf%>%
  mutate(Transect = as.factor(Transect))%>%
  group_by(Transect)%>%
  summarise()%>%
  st_cast("LINESTRING")

# getting study extent

study_ext <- st_bbox(points_sf)

ext <- st_as_sfc(study_ext)

# getting osm map

library(osmdata)

# Tamhini WLS

tamhini_raw <- opq(bbox = study_ext, timeout = 100)%>%
  add_osm_feature(key = 'boundary', value = 'protected_area')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
tamhini <- tamhini_raw$osm_multipolygons

st_crs(tamhini) <- "+init=epsg:4326"

# Roads

tamhini_area <- st_bbox(tamhini)

roads_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'highway')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
roads <- roads_raw$osm_lines

st_crs(roads) <- 4326

# Landuse

landuse_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'natural')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
forest <- landuse_raw$osm_polygons%>%
  filter(natural %in% c("wood", "tree"))

wood <- landuse_raw$osm_multipolygons%>%
  filter(natural == "wood")

water <- landuse_raw$osm_multipolygons%>%
  filter(natural == "water")

water_small <- landuse_raw$osm_polygons%>%
  filter(natural == "water")

st_crs(forest) <- 4326

# waterways

waterway_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'waterway')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
waterways <- waterway_raw$osm_lines

st_crs(waterways) <- 4326

# elevation 

library(elevatr)
library(rgdal)

tam_elev <- get_elev_raster(locations = tamhini, z = 8) 

tam_contours <- rasterToContour(tam_elev, levels = seq(0, 1000, 100))
# plotting site map

library(tmap)
library(cowplot) # for plot_grid() function - good to arrange multiple plots into a grid
library(grid)

# Study site

tm_shape(forest)+ #forested lands
  tm_polygons(col = "#81AF81", alpha = 0.5, border.col = NULL)+
tm_shape(wood)+ #forested lands
  tm_polygons(col = "#81AF81", alpha = 0.5, border.col = NULL)+
tm_shape(tam_contours)+
  tm_lines(col = "white", lty = "dotted")+
  tm_text(text = "level", auto.placement = T)+
tm_shape(water)+
  tm_polygons(col = "#77A1CB", border.col = NULL)+
tm_shape(waterways)+
  tm_lines(col = "#77A1CB", lwd = 2)+
tm_shape(roads)+
  tm_lines(alpha = 0.5, lwd = 1.5)+
  tm_text("ref", size = 0.3, remove.overlap = T)+
tm_shape(tamhini)+
  tm_borders(lwd = 2, col = "#165916", lty = "dashed")+
  tm_text("name", size = 0.5)+
tm_shape(lines, is.master = T, bbox = study_ext+c(-0.005, -0.005,0.005,0.005))+
  tm_lines(col = "red", lwd = 2)+
  tm_text("Transect", size = 0.5, auto.placement = T, remove.overlap = T)+
  tm_graticules(n.x = 5, lines = F, labels.inside.frame = T)+
  tm_scale_bar(position = c("left", "top"), text.size = 0.5, size = 0.5) + 
  tm_compass(position = c('right', 'top'), size = 2)+
  tm_layout(scale = 1.5, legend.position = c("left", "bottom"), bg.color = "#ffdb9c", asp = 1.5)
fig1a <- grid.grab()
# Tamhini WLS

tm_shape(forest)+ #forested lands
  tm_polygons(col = "#81AF81", alpha = 0.5, border.col = NULL)+
tm_shape(tam_contours)+
  tm_lines(col = "white", lty = "dotted")+
  tm_text(text = "level", auto.placement = T)+
tm_shape(water)+
  tm_polygons(col = "#77A1CB", border.col = NULL)+
  tm_text(text = "name", remove.overlap = F, size = 0.5)+  
tm_shape(water_small)+
  tm_polygons(col = "#77A1CB", border.col = NULL)+
  tm_text("name", remove.overlap = F)+  
tm_shape(waterways)+
  tm_lines(col = "#77A1CB", lwd = 2)+
tm_shape(roads)+
  tm_lines(alpha = 0.5, lwd = 1.5)+
tm_shape(tamhini, is.master = T, bbox = tamhini_area+c(-0.01, -0.01,0.01,0.01))+
  tm_borders(col = "#165916", lty = "dashed")+
  tm_text("name", size = 0.5)+
tm_shape(ext)+
  tm_borders(col = "red")+
  tm_scale_bar(position = c("left", "top"), text.size = 0.5) + 
  tm_layout(scale = 1, legend.position = c("left", "bottom"), bg.color = "#ffdb9c", asp = 1)
fig1b <- grid.grab()
# legend

tm_shape(forest)+
  tm_polygons(col = "#81AF81", alpha = 0.5, border.col = NULL)+
tm_shape(ext)+
  tm_borders(col = "red")+
  tm_add_legend(type = "fill",col = c("#81AF81", "#77A1CB", "#ffdb9c"), 
                labels = c("Forest", "Water", "Barren Land"), title = "Legend")+
  tm_add_legend(type = "line", col = c("black", "#165916", "#77A1CB", "red"), 
                labels = c("Roads", "WLS boundary", "Waterway", "Transect"), 
                lwd = c(3,3,3,3),
                lty = c(1, 2, 1, 1))+
  tm_layout(legend.only = T, asp = 1, legend.position = c(0.3,0.2), scale = 2)
legend <- grid::grid.grab()
# making inset map

tamext <- st_as_sfc(tamhini_area + c(-.5,-.5, .5, .5)) 

india <- read_sf("./Data/India-States.shp")

tmap_options(check.and.fix = TRUE)

tm_shape(india)+
  tm_polygons(col = "#ffdb9c")+
tm_shape(tamext)+
  tm_borders(lwd = 1, col =  "red")+
  tm_layout(scale = 1.5, bg.color = "#77A1CB", asp = 1)
fig1c <- grid.grab()
# making final map

library(gridExtra)

fig1 <- grid.arrange(fig1a, fig1b, fig1c, legend, layout_matrix = rbind(c(1, 1, 1),
                                                          c(1, 1, 1),
                                                          c(2, 3, 4)))

ggsave(fig1, filename = "./Figures and Tables/figure1.png", height = 8, width = 8)